knitr::opts_chunk$set(fig.width=5, fig.height=5) 

rm(list=ls())


options(scipen=999)

pkgs <- c("dplyr", "ggplot2", "stargazer", "car", "tidyverse", "extrafont", "haven", "coefplot", "plotrix", "sensemakr", "mice")
# A function to load all above packages. Install if they have not been installed.
usePackage <- function(p){
  for (pkg in p){
    if (!is.element(pkg, installed.packages()[,1]))
      install.packages(pkg, dep = TRUE, repos = "https://cloud.r-project.org/")
    require(pkg, character.only = TRUE)
  }
}
usePackage(pkgs)


d0 <- read.csv("Bram_CivicPulse_2023_05.csv")


d0$Republican <- 0
d0$Republican[d0$Demo_party_SC == "Republican" | d0$Demo_partylean_SC == "Republican Party"] <- 1
d0$Democrat <- 0
d0$Democrat[d0$Demo_party_SC == "Democrat" | d0$Demo_partylean_SC == "Democratic Party"] <- 1

d0$Independent <- as.numeric(ifelse(d0$Democrat == 1 | d0$Republican == 1, 0, 1))


d <- subset(d0, d0$Independent == 0 )

d$WardTreatment <- ifelse(d$Party1_EF != "officials", 1 ,0)
d$BoardSizeTreatment <- ifelse(d$Party2_EF != "officials", 1,0)

## FTs


d$FeelingDemCandidates <- case_when(
  d$Sentiment_GR_demcand == "No feeling at all" ~ .5,
  d$Sentiment_GR_demcand == "Somewhat favorable feeling" ~ .75,
  d$Sentiment_GR_demcand == "Very favorable feeling" ~ 1,
  d$Sentiment_GR_demcand == "Very unfavorable feeling" ~ 0,
  d$Sentiment_GR_demcand == "Somewhat unfavorable feeling" ~ .25,
  TRUE ~ NA_real_
)
d$FeelingDemCandidates <- as.numeric(d$FeelingDemCandidates)


d$FeelingDemCandidates <- case_when(
  d$Sentiment_GR_demcand == "No feeling at all" ~ .5,
  d$Sentiment_GR_demcand == "Somewhat favorable feeling" ~ .75,
  d$Sentiment_GR_demcand == "Very favorable feeling" ~ 1,
  d$Sentiment_GR_demcand == "Very unfavorable feeling" ~ 0,
  d$Sentiment_GR_demcand == "Somewhat unfavorable feeling" ~ .25,
  TRUE ~ NA_real_
)
d$FeelingDemCandidates <- as.numeric(d$FeelingDemCandidates)


d$FeelingDemVoters <- case_when(
  d$Sentiment_GR_demvoter == "No feeling at all" ~ .5,
  d$Sentiment_GR_demvoter == "Somewhat favorable feeling" ~ .75,
  d$Sentiment_GR_demvoter == "Very favorable feeling" ~ 1,
  d$Sentiment_GR_demvoter == "Very unfavorable feeling" ~ 0,
  d$Sentiment_GR_demvoter == "Somewhat unfavorable feeling" ~ .25,
  TRUE ~ NA_real_
)
d$FeelingDemVoters <- as.numeric(d$FeelingDemVoters)


d$FeelingRepCandidates <- case_when(
  d$Sentiment_GR_repcand == "No feeling at all" ~ .5,
  d$Sentiment_GR_repcand == "Somewhat favorable feeling" ~ .75,
  d$Sentiment_GR_repcand == "Very favorable feeling" ~ 1,
  d$Sentiment_GR_repcand == "Very unfavorable feeling" ~ 0,
  d$Sentiment_GR_repcand == "Somewhat unfavorable feeling" ~ .25,
  TRUE ~ NA_real_
)
d$FeelingRepCandidates <- as.numeric(d$FeelingRepCandidates)


d$FeelingRepVoters <- case_when(
  d$Sentiment_GR_repvoter == "No feeling at all" ~ .5,
  d$Sentiment_GR_repvoter == "Somewhat favorable feeling" ~ .75,
  d$Sentiment_GR_repvoter == "Very favorable feeling" ~ 1,
  d$Sentiment_GR_repvoter == "Very unfavorable feeling" ~ 0,
  d$Sentiment_GR_repvoter == "Somewhat unfavorable feeling" ~ .25,
  TRUE ~ NA_real_
)
d$FeelingRepVoters <- as.numeric(d$FeelingRepVoters)


d$FTOutPartyPeople <- ifelse(d$Republican == 1, d$FeelingDemVoters, d$FeelingRepVoters)
d$FTInPartyPeople <- ifelse(d$Republican == 1, d$FeelingRepVoters, d$FeelingDemVoters)

d$AffPolPeople <- d$FTInPartyPeople - d$FTOutPartyPeople

d$FTOutPartyPols <- ifelse(d$Republican == 1, d$FeelingDemCandidates, d$FeelingRepCandidates)
d$FTInPartyPols <- ifelse(d$Republican == 1, d$FeelingRepCandidates, d$FeelingDemCandidates)

d$AffPolPols <- d$FTInPartyPols - d$FTOutPartyPols

d$AffPol <- (d$AffPolPols + d$AffPolPeople)/2

## AFF POL GRAPHS
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)

# Convert sentiment responses to sentiment categories
d <- d %>%
  mutate(
    FeelingDemCandidates = case_when(
      Sentiment_GR_demcand == "Very unfavorable feeling" ~ "Very unfavorable",
      Sentiment_GR_demcand == "Somewhat unfavorable feeling" ~ "Somewhat unfavorable",
      Sentiment_GR_demcand == "No feeling at all" ~ "No feeling at all",
      Sentiment_GR_demcand == "Somewhat favorable feeling" ~ "Somewhat favorable",
      Sentiment_GR_demcand == "Very favorable feeling" ~ "Very favorable",
      TRUE ~ NA_character_
    ),
    FeelingDemVoters = case_when(
      Sentiment_GR_demvoter == "Very unfavorable feeling" ~ "Very unfavorable",
      Sentiment_GR_demvoter == "Somewhat unfavorable feeling" ~ "Somewhat unfavorable",
      Sentiment_GR_demvoter == "No feeling at all" ~ "No feeling at all",
      Sentiment_GR_demvoter == "Somewhat favorable feeling" ~ "Somewhat favorable",
      Sentiment_GR_demvoter == "Very favorable feeling" ~ "Very favorable",
      TRUE ~ NA_character_
    ),
    FeelingRepCandidates = case_when(
      Sentiment_GR_repcand == "Very unfavorable feeling" ~ "Very unfavorable",
      Sentiment_GR_repcand == "Somewhat unfavorable feeling" ~ "Somewhat unfavorable",
      Sentiment_GR_repcand == "No feeling at all" ~ "No feeling at all",
      Sentiment_GR_repcand == "Somewhat favorable feeling" ~ "Somewhat favorable",
      Sentiment_GR_repcand == "Very favorable feeling" ~ "Very favorable",
      TRUE ~ NA_character_
    ),
    FeelingRepVoters = case_when(
      Sentiment_GR_repvoter == "Very unfavorable feeling" ~ "Very unfavorable",
      Sentiment_GR_repvoter == "Somewhat unfavorable feeling" ~ "Somewhat unfavorable",
      Sentiment_GR_repvoter == "No feeling at all" ~ "No feeling at all",
      Sentiment_GR_repvoter == "Somewhat favorable feeling" ~ "Somewhat favorable",
      Sentiment_GR_repvoter == "Very favorable feeling" ~ "Very favorable",
      TRUE ~ NA_character_
    ),
    # Assign in-party and out-party sentiments
    InPartyPeople = ifelse(Republican == 1, FeelingRepVoters, FeelingDemVoters),
    OutPartyPeople = ifelse(Republican == 1, FeelingDemVoters, FeelingRepVoters),
    InPartyPols = ifelse(Republican == 1, FeelingRepCandidates, FeelingDemCandidates),
    OutPartyPols = ifelse(Republican == 1, FeelingDemCandidates, FeelingRepCandidates)
  )

# Combine data into long format
d_long <- d %>%
  select(InPartyPeople, OutPartyPeople, InPartyPols, OutPartyPols) %>%
  pivot_longer(cols = everything(), names_to = "Category", values_to = "Sentiment") %>%
  mutate(
    Group = case_when(
      str_detect(Category, "People") ~ "People",
      str_detect(Category, "Pols") ~ "Politicians"
    ),
    PartyType = case_when(
      str_detect(Category, "InParty") ~ "In-party",
      str_detect(Category, "OutParty") ~ "Out-party"
    )
  ) %>%
  filter(!is.na(Sentiment))

# Calculate counts and proportions
sentiment_counts <- d_long %>%
  group_by(Group, PartyType, Sentiment) %>%
  summarise(Count = n(), .groups = 'drop')

sentiment_totals <- sentiment_counts %>%
  group_by(Group, PartyType) %>%
  summarise(Total = sum(Count), .groups = 'drop')

sentiment_proportions <- sentiment_counts %>%
  left_join(sentiment_totals, by = c("Group", "PartyType")) %>%
  mutate(
    Proportion = Count / Total,
    SE = sqrt(Proportion * (1 - Proportion) / Total)
  )

# **Solution 1: Replace hyphens with underscores in PartyType**
sentiment_proportions <- sentiment_proportions %>%
  mutate(PartyType = str_replace(PartyType, "-", "_"))

# Pivot wider without hyphens in variable names
sentiment_diff <- sentiment_proportions %>%
  select(Group, Sentiment, PartyType, Proportion, SE) %>%
  pivot_wider(
    names_from = PartyType,
    values_from = c(Proportion, SE)
  ) %>%
  mutate(
    Difference = Proportion_In_party - Proportion_Out_party,
    SE_Diff = sqrt(SE_In_party^2 + SE_Out_party^2),
    CI_Lower = Difference - 1.96 * SE_Diff,
    CI_Upper = Difference + 1.96 * SE_Diff
  )

# Reorder Sentiment factor levels with new lines
sentiment_diff$Sentiment <- factor(sentiment_diff$Sentiment, 
                                   levels = c("Very unfavorable", "Somewhat unfavorable", "No feeling at all", "Somewhat favorable", "Very favorable"),
                                   labels = c("Very\nunfavorable", "Somewhat\nunfavorable", "No feeling\nat all", "Somewhat\nfavorable", "Very\nfavorable"))

# Increase font sizes
base_size <- 14

# Create the plot showing differences with error bars
plot <- ggplot(sentiment_diff, aes(x = Sentiment, y = Difference * 100, fill = Group)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), color = "black", alpha = 0.7) +
  geom_errorbar(aes(ymin = CI_Lower * 100, ymax = CI_Upper * 100), width = 0.2, position = position_dodge(width = 0.9)) +
  scale_fill_manual(name = "Group", values = c("People" = "#008080", "Politicians" = "#FF7F50")) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  theme_minimal(base_size = base_size) +
  theme(
    text = element_text(family = "LM Roman 10"),
    plot.title = element_text(hjust = 0.5),
    axis.text = element_text(size = base_size + 2),
    axis.title = element_text(size = base_size + 4),
    legend.title = element_text(size = base_size + 2),
    legend.text = element_text(size = base_size + 2),
    strip.text = element_text(face = "bold", size = base_size + 2),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
    legend.position = "bottom",
    strip.background = element_rect(fill = "gray90", color = "black", size = 0.5)
  ) +
  labs(
    x = "",
    y = "In-party rating minus out-party rating",
    title = "",
    fill = ""
  )

print(plot)
# Save the plot
ggsave(filename = "Ratings_of_Politicians_and_Public.jpeg", plot = plot, device = "jpeg", height = 6, width = 8, units = "in", dpi = 500)


## Demographics

# Create a function to extract the first year of the range and calculate age
age_from_range <- function(x) {
  if(grepl("earlier", x)) {
    return(2023 - 1920)
  } else if(grepl("later", x)) {
    return(2023 - 2006)
  } else {
    years <- as.numeric(unlist(strsplit(x, split = " - ")))
    return(2023 - years[1])
  }
}

# Apply the function to the age range column
d$Age <- sapply(d$Demo_age_DD, age_from_range)
d$Age <- ifelse(d$Age == 2122, NA, d$Age)

d$NonWhite <- as.numeric(ifelse(d$Demo_ethrace_MS_BIN == 1, 1, 0))

summary(as.factor(d$Demo_gender_SC))
d$Female <- ifelse(d$Demo_gender_SC == "Woman", 1, 0)

### -3 = Extremely Conservative; 3 = Extremely Liberal
summary(as.factor(d$Demo_ideo_SC))

d$Ideology <- case_when(
  d$Demo_ideo_SC == "Not sure" ~ .5,
  d$Demo_ideo_SC == "Moderate, middle of the road" ~ .5,
  d$Demo_ideo_SC == "Somewhat conservative" ~ .25,
  d$Demo_ideo_SC == "Somewhat liberal" ~ .75,
  d$Demo_ideo_SC == "Very liberal" ~ 1,
  d$Demo_ideo_SC == "Very conservative" ~ 0,
  TRUE ~ NA_real_
)
d$Ideology <- as.numeric(d$Ideology)


## Education: 1 = LT HS; 2 = High school; 3 = Trade school; 4 = Some college; 5 = college; 6 = some graduate school; 7 = grad degree
d$Education <- case_when(
  d$Demo_education_SC == "Less than high school" ~ 1,
  d$Demo_education_SC == "High school graduate" ~ 2,
  d$Demo_education_SC == "Technical/trade school" ~ 3,
  d$Demo_education_SC == "Some college" ~ 4,
  d$Demo_education_SC == "College graduate" ~ 5,
  d$Demo_education_SC == "Some graduate school" ~ 6,
  d$Demo_education_SC == "Graduate degree" ~ 7,
  TRUE ~ NA_real_
)
d$Education <- as.numeric(d$Education)


library(dplyr)
library(kableExtra)
# Calculate summary statistics with standard deviations and total N
demographic_summary <- d %>%
  summarise(
    Mean_Age = round(mean(Age, na.rm = TRUE), 2),
    SD_Age = round(sd(Age, na.rm = TRUE), 2),
    Proportion_NonWhite = round(mean(NonWhite, na.rm = TRUE), 2),
    Proportion_Female = round(mean(Female, na.rm = TRUE), 2),
    Mean_Education = round(mean(Education, na.rm = TRUE), 2),
    SD_Education = round(sd(Education, na.rm = TRUE), 2),
    Proportion_Democrat = round(mean(Democrat, na.rm = TRUE), 2),
    Proportion_Republican = round(mean(Republican, na.rm = TRUE), 2),
    Total_N = n()
  )

# Transpose the table
transposed_summary <- t(demographic_summary)
colnames(transposed_summary) <- c("Summary Statistics")

# Convert the transposed summary to a data frame and add row names
transposed_summary_df <- as.data.frame(transposed_summary)
transposed_summary_df$Variable <- c(
  "Mean Age",
  "SD Age",
  "Proportion NonWhite",
  "Proportion Female",
  "Mean Education",
  "SD Education",
  "Proportion Democrat",
  "Proportion Republican",
  "Total N"
)

# Generate LaTeX table
latex_table <- transposed_summary_df %>%
  select(Variable, `Summary Statistics`) %>%
  kable("latex", booktabs = TRUE, align = "c", col.names = c("Variable", "Summary Statistics")) %>%
  kable_styling(full_width = FALSE)

# Print the table
cat(latex_table)

## RESPONSE


d$WardResponse <- case_when(
  d$Ward_support_SC == "Strongly support" ~ 3,
  d$Ward_support_SC == "Support" ~ 2,
  d$Ward_support_SC == "Slightly support" ~ 1,
  d$Ward_support_SC == "Neither support nor oppose" ~ 0,
  d$Ward_support_SC == "Slightly oppose" ~ -1,
  d$Ward_support_SC == "Oppose" ~ -2,
  d$Ward_support_SC == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$WardResponse <- (d$WardResponse + 3) / 6
d$WardResponse <- as.numeric(d$WardResponse)

b1 <- lm(d$WardResponse ~ d$WardTreatment)
b1w <- lm(d$WardResponse ~ d$WardTreatment, weights = d$Weight)

###


d$SOWardDems <- case_when(
  d$Ward_dem_SC == "All" ~ 3,
  d$Ward_dem_SC == "Almost all" ~ 2,
  d$Ward_dem_SC == "Most" ~ 1,
  d$Ward_dem_SC == "About half" ~ 0,
  d$Ward_dem_SC == "Less than half" ~ -1,
  d$Ward_dem_SC == "Almost none" ~ -2,
  d$Ward_dem_SC == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOWardDems <- (d$SOWardDems + 3) / 6
d$SOWardDems <- as.numeric(d$SOWardDems)



d$SOWardReps <- case_when(
  d$Ward_rep_SC == "All" ~ 3,
  d$Ward_rep_SC == "Almost all" ~ 2,
  d$Ward_rep_SC == "Most" ~ 1,
  d$Ward_rep_SC == "About half" ~ 0,
  d$Ward_rep_SC == "Less than half" ~ -1,
  d$Ward_rep_SC == "Almost none" ~ -2,
  d$Ward_rep_SC == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOWardReps <- (d$SOWardReps + 3) / 6
d$SOWardReps <- as.numeric(d$SOWardReps)


d$SOWardInParty <- ifelse(d$Democrat == 1, d$SOWardDems, d$SOWardReps)
d$SOWardOutParty <- ifelse(d$Democrat == 1, d$SOWardReps, d$SOWardDems)

warddems <- lm(d$SOWardDems ~ d$Democrat )
wardreps <- lm(d$SOWardReps ~ d$Republican)


###


d$BoardSizeResponse <- case_when(
  d$Size_support_SC == "Strongly support" ~ 3,
  d$Size_support_SC == "Support" ~ 2,
  d$Size_support_SC == "Slightly support" ~ 1,
  d$Size_support_SC == "Neither support nor oppose" ~ 0,
  d$Size_support_SC == "Slightly oppose" ~ -1,
  d$Size_support_SC == "Oppose" ~ -2,
  d$Size_support_SC == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$BoardSizeResponse <- (d$BoardSizeResponse + 3) / 6
d$BoardSizeResponse <- as.numeric(d$BoardSizeResponse)

b2 <- lm(d$BoardSizeResponse ~ d$BoardSizeTreatment)
b2w <- lm(d$BoardSizeResponse ~ d$BoardSizeTreatment, weights = d$Weight)


###


d$SOBoardSizeDems <- case_when(
  d$Size_dem_SC == "All" ~ 3,
  d$Size_dem_SC == "Almost all" ~ 2,
  d$Size_dem_SC == "Most" ~ 1,
  d$Size_dem_SC == "About half" ~ 0,
  d$Size_dem_SC == "Less than half" ~ -1,
  d$Size_dem_SC == "Almost none" ~ -2,
  d$Size_dem_SC == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOBoardSizeDems <- (d$SOBoardSizeDems + 3) / 6
d$SOBoardSizeDems <- as.numeric(d$SOBoardSizeDems)



d$SOBoardSizeReps <- case_when(
  d$Size_rep_SC == "All" ~ 3,
  d$Size_rep_SC == "Almost all" ~ 2,
  d$Size_rep_SC == "Most" ~ 1,
  d$Size_rep_SC == "About half" ~ 0,
  d$Size_rep_SC == "Less than half" ~ -1,
  d$Size_rep_SC == "Almost none" ~ -2,
  d$Size_rep_SC == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOBoardSizeReps <- (d$SOBoardSizeReps + 3) / 6
d$SOBoardSizeReps <- as.numeric(d$SOBoardSizeReps)



d$SOBoardSizeInParty <- ifelse(d$Democrat == 1, d$SOBoardSizeDems, d$SOWardReps)
d$SOBoardSizeOutParty <- ifelse(d$Democrat == 1, d$SOWardReps, d$SOBoardSizeDems)

boardsizedems <- lm(d$SOBoardSizeDems ~ d$Democrat)
boardsizereps <- lm(d$SOBoardSizeReps ~ d$Republican)



## PERCENT SUPPORT:

# Calculate support proportion for Democrats
dem_support_proportionBS <- d %>% 
  filter(Democrat == 1 & BoardSizeTreatment == 1) %>%
  summarise(sum(BoardSizeResponse >= 0.8333333) / n())

rep_support_proportionBS <- d %>% 
  filter(Republican == 1 & BoardSizeTreatment == 1) %>%
  summarise(sum(BoardSizeResponse >= 0.8333333) / n())

dem_support_proportionW <- d %>% 
  filter(Democrat == 1 & WardTreatment == 1) %>%
  summarise(sum(WardResponse >= 0.8333333) / n())

dem_support_proportionW <- aggregate(WardResponse ~ Democrat + WardTreatment, data = d, function(x) sum(x >= 0.8333333) / length(x))

rep_support_proportionW <- d %>% 
  filter(Republican == 1 & WardTreatment == 1) %>%
  summarise(sum(WardResponse >= 0.8333333) / n())


## AVERAGE RESPONSE + EFFECT OF COMPETITION + IDEOLOGY + THREAT 

## CANNOT DO EFFECT OF COMPETITION BECAUSE OF RESTRICTED DATA

d$AverageResponse <- (d$WardResponse + d$BoardSizeResponse)/2
d$Treatment <- ifelse(d$Bram_assignment_EF == "treatment", 1, 0)
d$AvgInParty <- (d$SOBoardSizeInParty + d$SOWardInParty)/2
d$AvgOutParty <- (d$SOBoardSizeOutParty + d$SOWardOutParty)/2

d$PolarizedCounty0 <- (abs(d$County_2020pres_voteshare_dem_RESTRICTED - d$County_2020pres_voteshare_rep_RESTRICTED))
d$PolarizedCounty <- 1-d$PolarizedCounty0
d$Ideologue <- ifelse(d$Demo_ideo_SC == "Very liberal" | d$Demo_ideo_SC == "Very conservative", 1, 0)

#d$AboveMedianCompetition <- as.factor(ifelse(d$PolarizedCounty > median(d$PolarizedCounty, na.rm = TRUE), 1, 0))

model <- lm(d$AverageResponse ~ d$Treatment*d$PolarizedCounty)

compplot <- ggplot(d, aes(x = PolarizedCounty, y = AverageResponse, color = factor(Treatment))) +
  stat_smooth(method = "lm", se = TRUE, fullrange = TRUE) +
  scale_color_discrete(name = "Condition", labels = c("Control", "Treatment")) +
  scale_x_continuous(name = "County-level election competitiveness\n (0 = least, 1 = most competitive)", limits = c(0, 1), expand = c(0, 0)) +
  theme_minimal() +
  theme(
    text = element_text(family = "LM Roman 10"), 
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = .5)
  ) +
  labs(
    y = "Support for policies that reduce\n minority party competitiveness", 
    title = "Competitive elections and policy support"
  ) +
  geom_rug(sides = "b", aes(color = factor(Treatment)), position = "jitter")

# Display the plot
compplot


ggsave(filename="compplot.jpeg", plot=compplot, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)

### THREAT


d$AverageThreat <- (d$SOWardOutParty + d$SOBoardSizeOutParty)/2
#d$AverageInPartyThreat <- (d$SOWardInParty + d$SOBoardSizeInParty)/2


summary(lm(d$AverageResponse ~ d$AverageThreat*d$Treatment))

library(ggExtra)


compplotThreat <- ggplot(d, aes(x = AverageThreat, y = AverageResponse, color = factor(Treatment))) +
  geom_point(alpha = 0) +  # Invisible points
  stat_smooth(method = "lm", se = TRUE, fullrange = TRUE) +
  scale_color_discrete(name = "Condition", labels = c("Control", "Treatment")) +
  theme_minimal() +
  theme(
    text = element_text(family = "LM Roman 10"), 
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = .5)
  ) +
  labs(
    x = "Perceived threat from opponents\n (0 = least, 1 = most threatening)", 
    y = "Support for policies that reduce\n minority party competitiveness", 
    title = "Threat from opponents and policy support"
  )

# Add marginal density plots
compplotThreat_with_marginals <- ggMarginal(compplotThreat, type = "density", margins = "x", groupColour = TRUE, groupFill = TRUE)

# Display the plot
compplotThreat_with_marginals

ggsave(filename="compplotThreat.jpeg", plot=compplotThreat_with_marginals, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)


#d$Ideologue <- case_when(
#  d$Demo_ideo_SC == "Very liberal" ~ 1,
#  d$Demo_ideo_SC == "Very conservative" ~ 1,
#  TRUE ~ 0
#)
#d$Ideologue <- as.factor(d$Ideologue)



d$Treatment <- as.factor(d$Treatment)
d$Ideologue <- as.factor(d$Ideologue)


# Calculate means, standard error, and confidence intervals
#summary_stats <- d %>%
 # group_by(Ideologue, Treatment) %>%
  #summarise(mean_response = mean(AverageResponse, na.rm = TRUE),
   #         se = sd(AverageResponse, na.rm = TRUE) / sqrt(sum(!is.na(AverageResponse))),
    #        lower = mean_response - qt(0.975, sum(!is.na(AverageResponse))-1) * se,
     #       upper = mean_response + qt(0.975, sum(!is.na(AverageResponse))-1) * se)

summary_stats_manual <- data.frame()

# Unique combinations of Ideologue and Treatment
unique_combinations <- unique(d[, c("Ideologue", "Treatment")])

# Loop over each unique combination
for(i in 1:nrow(unique_combinations)) {
  
  # Subsetting the data
  subset_data <- subset(d, Ideologue == unique_combinations$Ideologue[i] & Treatment == unique_combinations$Treatment[i])
  
  # Calculate summary statistics
  mean_response <- mean(subset_data$AverageResponse, na.rm = TRUE)
  se <- sd(subset_data$AverageResponse, na.rm = TRUE) / sqrt(sum(!is.na(subset_data$AverageResponse)))
  lower <- mean_response - qt(0.975, sum(!is.na(subset_data$AverageResponse)) - 1) * se
  upper <- mean_response + qt(0.975, sum(!is.na(subset_data$AverageResponse)) - 1) * se
  
  # Populate the summary_stats_manual data frame
  summary_stats_manual <- rbind(summary_stats_manual, 
                                data.frame(Ideologue = unique_combinations$Ideologue[i],
                                           Treatment = unique_combinations$Treatment[i],
                                           mean_response = mean_response,
                                           se = se,
                                           lower = lower,
                                           upper = upper))
}


# Perform t-tests and tidy up the results
ttest_results <- d %>%
  group_by(Ideologue) %>%
  do(tidy(t.test(AverageResponse ~ Treatment, data = .))) %>%
  select(Ideologue, p.value)

# Create significance labels
ttest_results$sig_label <- case_when(
  ttest_results$p.value < 0.01 ~ "***",
  ttest_results$p.value < 0.05 ~ "**",
  ttest_results$p.value < 0.1 ~ "*",
  TRUE ~ "ns"
)


# Merge the t-test results back into the summary_stats data frame
summary_stats <- left_join(summary_stats_manual, ttest_results, by = "Ideologue")


interaction_test <- lm(AverageResponse ~ Ideologue * Treatment, data = d)
interaction_pvalue <- summary(interaction_test)$coefficients[4, 4]

# Assuming interaction_pvalue is the p-value for the interaction term
interaction_pvalue <- round(interaction_pvalue, 3)  # rounding for display purposes

# Determine the ymax for the significance bracket
ymax_bracket <- max(summary_stats$upper) + 0.1  # Adjust as needed

compplot_interact <- ggplot(summary_stats, aes(x = factor(Ideologue), y = mean_response, color = factor(Treatment), group = factor(Treatment))) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(width = 0.6)) +
  geom_point(position = position_dodge(width = 0.6), size = 4) +
  scale_x_discrete(labels = c("Weak ideology", "Strong ideologue")) +
  scale_color_discrete(name = "Condition", labels = c("Control", "Treatment")) +
  theme_minimal() +
  theme(
    text = element_text(family = "LM Roman 10"),
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = .5)
  ) +
  labs(
    x = "Ideological strength",
    y = "Support for policies that reduce\n minority party competitiveness",
    title = "Interaction between Ideology and Policy Support"
  ) #+   geom_hline(yintercept = c(0.345, 0.476), linetype = "dashed", color = "black")


# Display the plot
compplot_interact
ggsave(filename="compplot_interact.jpeg", plot=compplot_interact, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)




d$Ideologue <-as.numeric(d$Ideologue)

d$IdeologueAndCompetetive <- (d$Ideologue+d$PolarizedCounty)/2

model3 <- lm(d$AverageResponse ~ d$Treatment*d$IdeologueAndCompetetive)
compplot3 <- ggplot(d, aes(x = IdeologueAndCompetetive, y = AverageResponse, color = factor(Treatment))) +
  stat_smooth(method = "lm") +
  scale_color_discrete(name = "Condition", labels = c("Control", "Treatment")) +
  theme_minimal() +
  theme(text = element_text(family = "LM Roman 10"), 
        plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = .5)) +
  labs(x = "0 = Weak-ideologue / least competetive\n 1 = Strong ideologue / most competitive", 
       y = "Support for policies that reduce\n minority party competitiveness", 
       title = "Combination of ideology and local political competition")
compplot3
ggsave(filename="compplot3.jpeg", plot=compplot3, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)


compplotIntmodel<- lm(d$AverageResponse ~ d$Treatment*d$Ideologue)
compplotIntmodel2<- lm(d$AverageResponse ~ d$Treatment*d$PolarizedCounty)
compplotIntmodel3<- lm(d$AverageResponse ~ d$Treatment*d$AverageThreat)


### AFF POL

summary(lm(d$AverageResponse ~ d$AffPol*d$Treatment))

library(ggExtra)

compplotThreat <- ggplot(d, aes(x = AffPol, y = AverageResponse, color = factor(Treatment))) +
  geom_point(alpha = 0) +  # Invisible points
  stat_smooth(method = "lm", se = TRUE, fullrange = TRUE) +
  scale_color_discrete(name = "Condition", labels = c("Control", "Treatment")) +
  theme_minimal() +
  theme(
    text = element_text(family = "LM Roman 10"), 
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = .5)
  ) +
  labs(
    x = "Perceived affective polarization from opponents\n (0 = least, 1 = most polarized)", 
    y = "Support for policies that reduce\n minority party competitiveness", 
    title = "Affective polarization and policy support"
  )

# Add marginal density plots
affpolthreatplot <- ggMarginal(compplotThreat, type = "density", margins = "x", groupColour = TRUE, groupFill = TRUE)

# Display the plot
affpolthreatplot

ggsave(filename="affpolthreatplot.jpeg", plot=affpolthreatplot, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)




## MAKE THE DATA LONG 
model <- lm(d$AverageResponse ~ d$Treatment )


##Response_Type variables 
d$WResponse <- d$WardResponse
d$BResponse <- d$BoardSizeResponse

## Treatment Type variables
d$TreatW <- d$WardTreatment
d$TreatB <- d$BoardSizeTreatment



y <- dplyr::mutate(d, ID = row_number())

## DVs to common stem
y$DV_1 <- y$WResponse
y$DV_2 <- y$BResponse
y$DV_3 <- y$SOWardInParty
y$DV_4 <- y$SOWardOutParty
y$DV_5 <- y$SOBoardSizeInParty
y$DV_6 <- y$SOBoardSizeOutParty

## Reshape to long
longS <- reshape(y, direction = "long", 
                 varying = list(c("DV_1","DV_2","DV_3","DV_4","DV_5","DV_6")))
longS$DV <- longS$time
longS$response <- longS$DV_1

## Code treatment variable
longS$treat <- ifelse(longS$DV == 1, as.numeric(longS$TreatW), 
                      ifelse(longS$DV == 2, as.numeric(longS$TreatB), NA)
)
longS$treat <- longS$treat - 1
table(longS$treat)


## Code in vs out party
longS$inorout <- ifelse(longS$DV %in% c(3,5), 1, 
                        ifelse(longS$DV %in% c(4,6), 0, NA))

## Code ward or board size
longS$WorB <- ifelse(longS$DV %in% c(3,4), 1, 
                     ifelse(longS$DV %in% c(5,6), 0, NA))

## Model for experiment

#m_1 <- lm(response ~ treat + I(DV - 1), data = subset(longS, DV == 1 | DV == 2))
#summary(m_1)

#USE MODEL
summary(model)



## Model for SO beliefs
library(sandwich)
library(mice)

m_2 <- lm(response ~ inorout + WorB, data = subset(longS, DV %in% 3:6))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)

coefs_clustered <- cbind(coef(m_2), ses_2, t_2, (1 - pt(t_2, m_2$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered, 
          title = "SO Beliefs Model", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))

## MAKE MAIN PLOT

control_mean_1 <- coef(model)[1]
control_se_1 <- sqrt(diag(vcovCL(model)))[1]

treat_mean_1 <- coef(model)[1] + coef(model)[2]
treat_se_1 <- sqrt(control_se_1^2 + (vcovCL(model))[1, 2]^2)

# Calculate mean and standard errors for intercept and "inorout" variable in m_2 model
control_mean_2 <- coef(m_2)[1]
control_se_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))[1]

inorout_mean_2 <- coef(m_2)[1] + coef(m_2)[2]
inorout_se_2 <- sqrt(control_se_2^2 + (vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6]))[1, 2]^2)

# Create table
table_data <- data.frame(
  variable = c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support"),
  mean = c(control_mean_1, treat_mean_1, inorout_mean_2, control_mean_2),
  se = c(control_se_1, treat_se_1, inorout_se_2, control_se_2)
)

mean_data3 <- as.data.frame(table_data)
mean_data3$variable <- factor(mean_data3$variable, levels = c("Control","Support the proposal", "Beliefs about in-party support",  "Beliefs about out-party support"))

mean_data4 <- head(mean_data3, nrow(mean_data3)-2)


# Create the bar plot with error bars
soplotcp <- ggplot(mean_data3, aes(x = variable, y = mean, fill = variable)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), aes(fill = variable)) +
  scale_fill_manual(values = c("#5aa9e6", "#74c67a", "#ff9866", "#8e76c3")) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black", position = "identity") +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "Local policymaker support for policy proposals\n that would reduce the power of the party in the minority",
       x = "", y = "Mean Response") + 
  scale_x_discrete(labels = c("Control:\n Non-partisan Support", "Treatment:\n Partisan Support", "Perceived\n In-Party Support", "Perceived\n Out-Party Support")) +  
  theme( axis.line.y = element_blank(), axis.ticks.y=element_blank())+   
  theme(text = element_text(size=10, family="LM Roman 10")) +  
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5, family = "LM Roman 10"))  +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  theme(legend.position = "none") +
  scale_y_continuous( labels = scales::percent_format(accuracy = 1),  limits =  c(0, 1))

soplotcp

ggsave(filename="soplotcp.jpeg", plot=soplotcp, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)


# Create the bar plot with error bars
soplotExpcp <- ggplot(mean_data4, aes(x = variable, y = mean, fill = variable)) +
  geom_point(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "Local policymaker support for policy proposals\n that would reduce the power of the party in the minority",
       x = "",
       y = "Mean level of support (50% = neither support nor oppose)") +
  scale_x_discrete(labels = c("Control -\n Support for proposal\n (No partisan information)", 
                              "Treatment -\n Support for proposal\n (Benefits in-party)")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
  theme_bw() +
  theme(axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        text = element_text(size=10, family="LM Roman 10"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "none")

soplotExpcp

ggsave(filename="soplotExpcp.jpeg", plot=soplotExpcp, device="jpeg", height=5, width=5, units="in", dpi=500)

##



## PARTISAN HETEROGENEITY NOT USED

if (!require(broom)) {
  install.packages('broom')
}


d$AffectivePolarization <- d$AffPol
# Creating a linear model
model1 <- lm(AverageResponse ~ Treatment + Republican + AffectivePolarization, data = d)

# Summarize the model results into a data frame
df_95 <- broom::tidy(model1, conf.int = TRUE, conf.level = 0.95)
df_90 <- broom::tidy(model1, conf.int = TRUE, conf.level = 0.90)

# Merge the two data frames
df <- merge(df_95, df_90, by = c("term", "estimate", "std.error", "statistic", "p.value"), 
            suffixes = c(".95", ".90"))

# Remove the intercept
df <- df[df$term != "(Intercept)",]

# Plotting
cplot1cp <- ggplot(df, aes(y = term, x = estimate)) +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(color = "black") +
  geom_segment(aes(x = conf.low.95, xend = conf.high.95, y = term, yend = term), size = 0.5, color = "black") +
  geom_segment(aes(x = conf.low.90, xend = conf.high.90, y = term, yend = term), size = 1, color = "black") +
  labs(title = "Percent increase in\n support for policies\n subverting minority power", 
       x = "", 
       y = "") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 2), limits = c(-.1, .3)) +
  theme_bw() +
  theme(text = element_text(size = 10, family = "LM Roman 10"),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position = "none")
cplot1cp

ggsave(filename="coefplotcp.jpeg", plot=cplot1cp, device="jpeg", height=5, width=5, units="in", dpi=500)

longS$outorin <-  ifelse(longS$inorout == 0, 1, 0)

m_2 <- lm(response ~ outorin + Republican + AffPol + WorB, data = subset(longS, DV %in% 3:6))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)

coefs_clustered_m2 <- cbind(coef(m_2), ses_2, t_2, (1 - pt(t_2, m_2$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered_m2) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered_m2, 
          title = "Model 2", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))


## GRAPH:
# Ensure the necessary packages are installed
if (!require(broom)) {
  install.packages('broom')
}
if (!require(multiwayvcov)) {
  install.packages('multiwayvcov')
}
if (!require(dplyr)) {
  install.packages('dplyr')
}
if (!require(scales)) {
  install.packages('scales')
}

# Load the packages
library(broom)
library(multiwayvcov)
library(dplyr)
library(scales)

# Subset the data
sub_data <- subset(longS, DV %in% 3:6)

# Create a linear model
m_2 <- lm(response ~ outorin + Republican + AffPol + WorB, data = sub_data)

# Get the variance-covariance matrix with clustered standard errors
vcov_clustered <- cluster.vcov(m_2, sub_data$id)

# Get standard errors
clustered_se <- sqrt(diag(vcov_clustered))

# Create a data frame of results
df <- tidy(m_2)
df$clustered_se <- clustered_se

# Compute confidence intervals
df$conf.low.95 <- df$estimate - 1.96 * df$clustered_se
df$conf.high.95 <- df$estimate + 1.96 * df$clustered_se
df$conf.low.90 <- df$estimate - 1.645 * df$clustered_se
df$conf.high.90 <- df$estimate + 1.645 * df$clustered_se

# Remove the intercept and WorB from the data frame
df <- df %>% filter(term != "(Intercept)", term != "WorB")

# Relabel the terms
df <- df %>%
  mutate(term = case_when(
    term == "outorin" ~ "Out-party vs. in-party",
    term == "AffPol" ~ "Affective Polarization",
    TRUE ~ term
  ))

# Plotting
cplot1cpso <- ggplot(df, aes(y = term, x = estimate)) +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(color = "black") +
  geom_segment(aes(x = conf.low.95, xend = conf.high.95, y = term, yend = term), size = 0.5, color = "black") +
  geom_segment(aes(x = conf.low.90, xend = conf.high.90, y = term, yend = term), size = 1, color = "black") +
  labs(title = "Predictors of beliefs about who\n would support policies\n that subvert minority power", 
       x = "", 
       y = "") +
  scale_x_continuous(labels = percent_format(accuracy = 1), limits = c(-.1, .3)) +
  theme_bw() +
  theme(text = element_text(size = 10, family = "LM Roman 10"),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position = "none")
cplot1cpso
ggsave(filename="coefplotcpso.jpeg", plot=cplot1cpso, device="jpeg", height=5, width=5, units="in", dpi=500)

### DEMOGRAPHICS MODEL


modeldem <- lm(AverageResponse ~ Treatment + Republican + AffPol + NonWhite + Education  + Female + Age , data = d)
summary(modeldem)

### partisan heterogeneity model for table 7

model <- lm(d$AverageResponse ~ d$Treatment*d$Republican)


